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ABSTRACT 

The  unified  theory  developed  by  Henrici  [3]  for  one-step  methods  is 
extended  to  the  more  general  case  where  the  order  of  the  system  of  difference 
equations  can  exceed  the  order  of  the  system  of  ordinary  differential  equa- 
tions.  The  analysis  is  applicable  to  every  fixed-stepsize  fixed-formula 
method  known  to  the  author.   For  many  of  these  methods  the  concept  of  consis- 
tency is  inadequate.   A  more  appropriate  concept,  termed  quasi-consistency, 
is  introduced. 


1.   Introduction .  Most  numerical  methods  for  solving  initial  value 
problems  in  ordinary  differential  equations  discretize  the  system  of  differential 
equations  and  then  solve  the  resulting  system  of  difference  equations.   Henrici 
[3]  has  developed  the  theory  for  one-step  methods.   This  paper  studies  the  more 
general  case  where  the  order  of  the  system  of  difference  equations  may  exceed 
the  order  of  the  system  of  differential  equations.   Just  as  ordinary  differential 
equations  are  more  easily  studied  when  written  as  a  first  order  system,  it  is 
likewise  profitable  to  study  difference  equations  as  a  first  order  system: 

u  =  S  u    +  h  ±(t     ,  ,  u   ;  h). 

-n     -n-1     ■*■  n-1  -n-1 

This  one-step  formulation  is  quite  powerful,  for  it  renders  the  analysis  more 
transparent . 

Two  basic  properties  of  methods  are  of  practical  concern:   convergence, 
which  means  that  the  accumulated  discretization  error  goes  to  zero  as  the  step- 
size  h  ■*■  0,  and  stability,  which  ensures  that  the  accumulated  roundoff  error 
does  not  grow  "too  fast"  as  h  •+   0.   Some  authors  (for  example,  Chartres  and 
Stepleman  [l])  combine  the  concepts  of  convergence  and  stability.   This  approach 
is  cumbersome,  and  it  does  not  reflect  current  practice  in  the  sense  that 
generally  the  roundoff  error  does  not  tend  to  zero  like  some  positive  power  of 
h. 

Convergent  one-value  (one-step)  methods  generate  difference  equations 
which  are  discrete  analogs  of  the  differential  equation,  and  as  a  result  the 
stability  properties  of  the  differential  equation  are  preserved  by  the  difference 
equation  for  small  enough  h.   For  multivalue  methods  the  difference  equations 
are  not  analogs  of  the  differential  equations  and  so  the  stability  properties 
of  the  differential  equation  might  not  be  faithfully  represented.   In  order  to 
ensure  that  this  is  not  the  case,  we  ought  to  require  that  a  multivalue  method 
be  "as  stable  as"  a  one-value  method.   The  notion  of  being  "as  stable  as"  is 


defined  in  §3  in  a  quite  natural  way.   And  then  in  §U  it  is  demonstrated  that 
being  as  stable  as  a  one-value  method  is  equivalent  to  requiring  that  the  matrix 
S  have  all  of  its  eigenvalues  inside  the  unit  circle  except  for  an  eigenvalue 
at  1  whose  index  is  one  and  whose  multiplicity  is  equal  to  the  number  of  initial 
values  required  for  the  differential  equation. 

There  are  common  methods,  like  Milne's  method,  which  are  not  as  stable 
as  a  one-step  method  but  which  can  be  expressed  in  a  form  for  which  they  satisfy 
the  strict  root  condition:   all  the  roots  of  the  minimal  polynomial  of  S  are 
inside  the  unit  circle  except  possibly  for  a  root  at  1 .   It  is  determined  that 
for  methods  satisfying  the  strict  root  condition  convergence  is  equivalent  to 
a  property  termed  quasi-consistency.   If  d  ,  n  =  0(l)N,  denotes  the  local  dis- 
cretization error,  then  quasi-consistency  means  that 

max    i  ,   i    _  , 

_  d    -*  0    as    h  ■*■  0 

0  £  n  £  N  '  — n  ' 

and 

max    j  E(    +cL+"'+d)|-»-0    as    h  +  0 

0  $  n  £  N  '   -0   -^         -n   ' 

where  E  =  S    .      On  the  other  hand  consistency  usually  means  that 

l<^l    +    |d1'    +    ,*'   +    |cLl~>0  as  h  +  0 

(see  Chartres  and  Stepleman  [l],  Stetter  [7»  pp.  55  75]). 

It  is  shown  that  the  order  of  convergence  is  always  equal  to  the  order 
of  quasi-consistency.   For  an  often-studied  class  of  methods  which  Stetter  [7, 
p.  310]  calls  straight  multi step  methods ,  quasi-consistency  is  equivalent  to 
consistency.   This  is  not  true  for  more  general  methods.   For  example,  Nordsieck'i 
[5]  method  is  quasi-consistent,  and  hence  convergent,  of  order  6  and  yet  only 
consistent  of  order  5-   Other  examples  are  given  in  §5.   To  have  convergence 
of  order  p,  it  suffices  to  have  that  d  =  0(h  )  and  that  the  essential  local 


— n 


truncation  error  E  d  =  0(h   ).   In  the  case  of  linear  multistep  methods  E  d 

is  just  the  local  truncation  error  divided  by  8„  +  8  +  •"  +  B   (cf.  Gear  [2, 

U    1  q 

p.  118]). 

In  §6  the  leading  term  in  the  asymptotic  expansion  of  the  global  error 
of  a  method  satisfying  the  strict  root  condition  is  determined  in  terms  of  the 
leading  term(s)  in  the  asymptotic  expansion  of  the  local  discretization  error. 


2.   Fixed-stepsize  methods.   Let  s  be  a  fixed  positive  integer.   We 
consider  a  fixed  class  of  systems  of  s  first  order  differential  equations 
(2.1a)  y'  =  f(t,  y) 

where  f  is  continuous  in  t  and  uniformly  Lipschitz  continuous  in  y  for 

s 

0  S  t  <  1,  y€=R.   For  any  initial  condition 

(2.1b)  y(0)  =  yQ 

a  unique  solution  y(t)  is  guaranteed.  As  an  example,  the  class  of  systems 

under  consideration  might  consist  of  those  equations  having  the  form 

(y1)'  =  y2  , 

(y2)'  =  f(t,  y\  y2)   . 

Let  q  and  k  be  fixed  positive  integers  with  k  >,   s.   The  methods  under 
consideration  have  three  components: 

(i)    a  mapping  which  to  each  f  assigns  a  starting  procedure 

o(y;  h)^-Rk 
continuous  for  y€^  RS  and  0  £  h  £  h  where  0  <  h  <,   1/q. 
( ii )   a  mapping,  called  a  formula,  which  to  each  f  assigns  a  forward 
step  procedure 

S  u  +  h  _^(t,  u;  h)  G  Rk 
where  S  is  a  k  x  k  matrix  independent  of  f  and  the  increment 
function  _£  is  uniformly  Lipschitz  continuous  in  u  and  continuous 
in  t  and  h  for  0$  t  Jl,  ugH  ,  and  0  S  h  S  h  . 
(iii)   a  mapping  which  to  each  solution  y(t)  of  a  problem  assigns  a 
correct  value  function 

v_(t;  h)  =  A  y(t)  +  h  z.(t;  h)  £  Rk 
where  A  is  a  k  x  s  matrix  of  rank  s  independent  of  y(t)  and 
_z(t;  h)  is  bounded  for  0  $  (q-l)h  S  t  S  1. 


Stetter  [7,  p.  7]  introduces  a  similar  mapping,  which  he  denotes  by  A  . 

n 


Note  that  y(t)  can  be  recovered  "by  means  of 

y(t)  =  (AT  A)-1  AT  v_(t;  0). 

Given  a  problem  and  an  h  ^  h  ,  a  method  generates  a  uniform  grid  and 

a  discrete  problem  on  that  grid.   The  grid  is 

(q  -  l)h  =  tQ  <  t1  <  ...  <  tH  S  1 

where  t  =  t   +  nh  and  N  is  the  largest  integer  such  that  t   £  1.   Since  there 

is  no  danger  of  confusion  we  do  not  exhibit  the  dependence  on  h  of  such 

quantities  as  t  and  N.   The  discrete  problem  on  the  grid  is 

(2.2a)        Uq  =  a(yQ;  h), 

(2.2b)        u  =  S  u   .  +  h  j)(t   .  ,  u  .  ;  h)s    n  =  l(l)N, 
— n     — n-1     ■*■  n-1  — n-1 

which  can  be  solved  to  yield  a  numerical  solution 

U=  (uq,  u1}  ...,  %)T. 

The  discrete  problem  (2.2)  is  normalized  so  that  the  unknown  data  vector  u 

— n 

appears  on  the  left  hand  side.   This  normalization  serves  to  define 
o_,  S,  and  i£_  uniquely. 

The  correct  value  function  y_(t;  h)  assigns  a  meaning  to  the  numerical 
solution  U  by  means  of 

1=  (zo>  zv  •••>  i^)T 

where  y  =  y_(t  ;  h) .   Hence  the  global  discretization  error  is  defined  to  be 


U  -  Y . 

Remark.      An  unfortunate  aspect  of  methods  with  k  >  s  is  that  the  data 

within  u  are  normally  inconsistent  with  the  differential  equation.   For 

T 
example,  if  u  =  [u  ,  u   ]  ,  then  there  is  usually  no  solution  y(t)  of  the 

differential  equation  such  that  y(t  )=  u  and  y(t   .,  )  =  u 

n    n        n-1     n-1 

Example    I.      A  straight  q-step  method  (Stetter's  [7,  p.  310]  terminology) 

determines  one  new  value  for  each  step: 

u  =  S  u    +  h  en  i|i(t  _ ,  u  ,  ;  h) 
-n     -n-1     —1    n-1  —n-1 


where 


S  = 


u      = 
— n 

[V  Vi'  • 

•  >   u 
n 

-q+lJ    ' 

"ai 

-a2       ... 

-a 

q- 

i       q 

1 

0 

0 

0 

0 

l  k 

0 

0 

and 

e±   =  [1,  0,  ...,  0]T. 

T 
Also  v_(t;  h)  =  [y(t),  y(t-h),  ...,  y(t  -  (q-l)  h) ]  .   For  the  special  case  of 

a  linear  q-step  method 


i|>(t,  u;  h)  =  gQ  4>  +  i     0  f(t-j  h,  uJ) 


where  <j>  solves  the  equation 

4>  =  f(t,  h  U  +  {  -In(-a.  uJ  +  h  6.  f(t-j  h,  uj))}). 
Example  2.      Perhaps  the  most  important  class  of  methods  which  are  not 
straight  multistep  methods  are  P(EC)  multistep  methods.   For  example  a  PEC 
method  with  a  second  order  Adams-Bashforth  predictor  and  a  third  order  Adams- 
Moulton  corrector  can  be  written 


u 


u* 


n-1 


10  0 
10  0 
0     1      0 


u 


n-1 


n-1 


u 


n-2 


+  h 


5/12     2/3     -1/12 
0  3/2      -1/2 


■n-1 


"n-2 


where 


,P 
n-2 


f-   „   =   f(t_    oJ   u?  J, 


n-29     n-2' 

f9  .    =   f(t     .  ,   uP  ,  ), 
n-1  n-1       n-1 

>P 


=  f(t    ,  u     ,    +  |h  fP  ,    -  i-h  fP  _). 
n  n'      n-1        2  n-1        2  n-2 


T 


Naturally  y_(t;   h)   =    [y(t),   y(t),   y(t-h)]    .      Alternatively  it   could  be  written 


u 
n 

1 

2/3 

-1/12 

n-1 

h  u' 
n 

= 

0 

0 

0 

h  u1  . 
n-1 

h  u1  . 
n-1 

0 

1 

0 

h  u' 
n-2 

5/12 


•  h  f(t    ,  u    ,    +  |  h  u'       -  \h  u'     ) 

n       n-1        2         n-1       2         n-2 

with  y_(t;   h)   =   [y(t),   h  y'(t),   h  y'(t-h)]T. 

Example   3.      Butcher's  three-stage   formula  of  order   four    (cf.    Stetter 
[7,    p.    275])    is   given  by 


M              — 

*" 

u 
n 

0 

1 

Un-1 

+  h 

h  Un-l/2 

0 

0 

h  Un-3/2 

1/6     2/3     1/6 
0  10 


"n-1 


n-1/2 

f5 

n 


where 


rn-l        "   f(tn-l-   "„-!>' 


n-1/2 


=   f(t 


u 


+  fh   f 


=.    n  t 


n-1/2'      n-1        k  n-1  "   h     n-3/2 


), 


=   f(t    ,   u         +  2h  f     .  ._   -  2h   f     _    +  h  u'    _/0] 
n       n-1  n-1/2  n-1  n-3/2 


and  by  y_(t;   h)   =    [y(t),   h  y'(t-h/2)]     . 

Example   4.      A  modified  linear  multistep  method  equivalent   to  a  fourth 
order  Adams -Moult  on  method   (Gear    [2,   p.    150]   with  M  =  <=°)    is   given  by 


u 
n 

1/2 

1 

1/2 

1/8 

u  . 
n-1 

3/8 

h  u' 
n 

0 

0 

0 

0 

h  u'  , 
n-1 

1 

n-1 

1/2 

1/3 

1/2 

5/2U 

Va 

+  h 

-1/2U 

h  u'  , 
n-1 

0 

1 

0 

0 

hUn-2 

0 

ij^(t     ,  ,   u     n  ;   h) 
n-1     -n-1 


where  y_(t;h)  =  [y(t),  h  y*(t),  y(t-h),  h  y'(t-h)]T  and  i|;(t,  u.i  h)  =  «p  the 


solution  of 


ip  =  f(t,  \  u1  +  u2  +  |  u3  +  |  uU  +  |  h  ^) 


i     i  k 

3.      Convergence  and  stability.      Let    | * |    be  a  norm   for  R    ,   and  define 

I  I    tt       v    I  I  max  I  I 

II  H-I  IL-  0«  n  *  N   I   Hfc-Zn   I- 

DEFINITION  3.1.   A  method  is  said  to  be  convergent   for  the  problem 
(2.1)  if 

II  U-Y_llOo")'0ash">'0- 
If  the  method  is  convergent  for  all  infinitely  smooth  problems  (2.1),  then  it 
is  simply  said  to  be  convergent. 

In  practice  U  cannot  be  computed  exactly;   we  must  take  into  account 
the  local  computing  error  which  arises  from  doing  operations  in  finite 
precision.   Hence  we  consider  the  perturbed  numerical  solution 
LL=  (Hq>  Has  •••>  Ujj)   defined  by 
Uq.  =  ^0;  h)  +  £o' 


u 
— n 


=Su  .  +  h  *(t  ,  ,  u  _ ;  h)  +  r  ,  n  =  l(l)N, 
-n-1     *-  n-1  -n-1       — n 

m  a 

for  any  perturbation  R_  =  (r^,  r_  ,  .  .  .,  r    )    .   The  dependence  of  U  on  R  is  not 
made  explicit  in  our  notation. 

Various  definitions  of  stability  as  an  absolute  concept  have  been 
proposed  but  most  of  them  seem  somewhat  arbitrary.   (An  exception  is  the  use 
of  Spijker's  [6]  norm  to  define  stability  for  straight  multistep  methods.)   It 
is  more  natural  to  define  stability  as  a  relative  concept. 

DEFINITION  3.2.  Method  2  is  said  to  be  at   least  as  stable   as  method  1 
if  for  each  problem  there  is  a  fixed  number  K  such  that 

II  "(2)  -u(2)  il  *  k  ||  i(1)  -u(1)  il 

for  any  h  $  h  and  any  perturbation  R  ^  IT  +    .   If ,  in  addition,  method  2 
is  at  least  as  stable  as  method  1,  then  the  methods  are  said  to  be  equally 
stable ;  otherwise  method  1  is  said  to  be  move  stable   than  method  2. 


10 

Since  good  stability  properties  are  important  for  practical  algorithms, 
it  is  desirable  that  a  method  satisfy  the  most  stringent  stability  requirement 
that  is  compatible  with  convergence. 

DEFINITION  3.3.  A  convergent  method  is  optimally  stable   if  there  is 
no  convergent  method  which  is  more  stable. 

Note.      In  general  it  is  not  possible  to  compare  the  stability  of  two 
arbitrary  methods,  and  hence  two  optimally  stable  methods  are  not  necessarily 
equally  stable. 

Remark.      Spijker's  [6]  definition  of  optimal  stability  differs  from 
ours  in  that  it  requires  an  optimally  stable  method  to  be  at  least  as  stable 
as  every  other  method  in  the  class  of  methods  under  consideration. 

It  is  shown  in  §U  that  the  optimally  stable  methods  are  those 
convergent  methods  that  satisfy  the  very  strict  root  condition,  which  is 
defined  below. 

DEFINITION  3.^.  A  method  is  said  to  satisfy  the  root  condition    (RC) 
if  the  roots  of  the  minimal  polynomial  of  S  are  either  inside  the  unit  circle 
or  on  the  unit  circle  and  simple.  A  method  satisfies  the  strict  root  condition 
(SRC)  if  the  roots  of  the  minimal  polynomial  are  inside  the  unit  circle  except 
possibly  for  a  simple  root  at  1.   A  method  satisfies  the  very  strict  root 
condition    (VSRC)  if  it  satisfies  the  SRC  and  if  S  has  at  most  s  eigenvalues 
equal  to  1. 

It  is  shown  by  Theorem  3.9  that  the  matrix  S  must  have  at  least  s 
eigenvalues  equal  to  1  in  order  for  the  method  to  be  convergent;  however  any 
additional  eigenvalues  on  the  unit  circle  impart  a  degradation  to  the  stability 
of  the  method. 

There  is  no  increase  in  the  complexity  of  the  theory  if  instead  of 


11 


assuming  the  VSRC  we  merely  require  that  the  strict  root  condition  be 
satisfied.   Furthermore,  it  is  then  possible  to  treat  marginally  useful  methods 
like  Milne's  method.  Hence  the  theory  is  developed  for  methods  satisfying  the 
SRC. 

The  following  lemma  gives  the  best  possible  qualitative  bound  on  the 
difference  U  -  U  in  terms  of  R  for  methods  satisfying  the  RC.   First  we 
introduce  some  notation.   For  any  k  x  k  matrix  T  let  [T]  denote  the  (N+l)  x 
(N+l)k  matrix 

I 
-T  I 


-T  I 


Then  the   quantity 


max 


J        -Moo        0  :?  n   $  N    '    j=0 


n       n-1 
L,  T  r 


■J 


(N+l)k 


is  a  norm  for  R 

LEMMA  3.5-     Given  a  method  satisfying  the  RC  and  a  problem   (2.1),   there 
exist  positive  constants  c  and  C  such  that  for  any  h  $  h„  and  any  R  €R 

S    M    U  -  U   I  I      *  C 

O  I    I       —  I     I   00 


,(N+l)k 


c    II    [S]"1  R 


[S]_1  R 


Proof.      Set    6     =  u     -  u   .      Then 
-n       -n       -n 


6     =SS         +h6  +r 

— n  —n-1  —n-1       — n 


(3.1) 
where 

5     n=lk.("t     n>u     -,;h)-\J;(t     ,,u     ,;h) 
—n-1       *-    n-1     -n-1  *-    n-1     -n-1 

By  assumption  there  exists  a  constant  L  such  that 


4-i  I  s  L 


4-i 
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Solving  the  difference  equation  (3.1)  gives 
(3.2) 


n   n-1    ~       n   n-1 


By  the  root  condition  there  exists  a  constant  B  such  that 

i  «n 


(3.3) 

Thus  (3.2)  "becomes 


L   <:   B. 


n-1 
6  $   h  B    .En        6, 

-n  '  j=o   '  -j 


[S]"1   R 


By  induction  on  n  it   follows  that 


6 
— n 


-1 


$    (1   +  h  B)n    ||    [S]"1  R    ||    .<   e*    ||    [S]    1  R    |  |, 


which  proves  the  second  inequality  in  the  lemma.   From  (3.2)  and  (3.3) 


,n-J 


max 


i  j£o  S—  r.  |  S  (1  +  n  h  B)  Q  ^  „  |  ^  |;  «  (1  +  B)  |  |  U  -  U  ||„, 

and  thus   I |  S  _1  R  j  !   $  (l  +  B)  | I  U  -  U  j   .   Q.E.D. 

II  —  I  ico  'II  —     —  I  loo 

When  a  norm,  such  as-  |  I  [S]    •  M  ,  permits  an  upper  hound  on  U  -  U, 
then  the  method  is  said  to  he  stable  with  respect  to   that  norm,  and  if  the 
norm  also  permits  a  lower  bound  on  U  -  U,  then  the  norm  is  said  to  he  a 
minimal  stability  functional   for  the  method.   These  ideas  are  discussed  in 
more  detail  by  Spijker  [6]  and  Stetter  [7,  pp.  81-8U] . 

Of  considerable  importance  to  the  theory  is  the  component  of  S  (cf. 
Lancaster  [U,  p.  162])  corresponding  to  the  eigenvalue  1,  which  we  denote  by 
E.   Assuming  S  satisfies  the  RC,  we  have  that 

E  =  p"(l)_1  p"(S) 
where  p(£)  =  (g  -  l)p(^)  is  the  minimal  polynomial  of  S.   If  S  does  not  have 
an  eigenvalue  at  1,  then  we  define  E  =  0.   It  is  easily  verified  that  E  is 
idempotent  and  that  S  E  =  E  =  E  S  whence 

Sn  =  E  +  (S  -  E)n,     n  =  1,  2,  ...   . 
If  S  satisfies  the  SRC,  the  eigenvalues  of  S  -  E  are  inside  the  unit  circle. 
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The  following  theorem  is  for  methods  satisfying  the  SRC,  and  it  gives 

bounds  on  the  difference  U  -  U  in  terms  of  the  more  useful  norm 

I  I  r-,-1-1  „  I  I   _    max    I    n~1         I 

LE J   R     =  „       .T   E  .Z.  r.  +  r 
l|LJ   -i'oo   0  s  n  £  N  '    j=0  -j   -n  ' 

THEOREM  3.6.  Given  a  method  satisfying  the  SRC  and  a  problem   (2.1), 

there  exist  positive  constants   c'  and   C  such  that  for  any   h  £  h  and  any 
EeR(H+l)k 

c1     [E]_1  R  | !   £   |  U  -  U     £  C     [E]-1  R 
Proof.      The  theorem  follows  from  Lemma  3.5  if  it  can  be  shown  that 
[El   R  I    can  be  bounded  by  some  constant  times  '   [S]   R     and  vice 
versa.   First  of  all 

II  [E]"1  R  II.-  ||  [E]-1  [S]  [S]"1  R  || 
=  ||  [S  -  E]  [S]_1  R  II. 

=   0STSHl    Jo<tS-E])n.    ([Sl^R).    | 

N  i 

«0*T«N     jfio    I    (SS-^'nJ    I     II    tSl^R    IL 

=    (1   +    I    S   -   E    | )  [S]"1  R 

I  III  I  1       QQ 

and  secondly 

II    [S]"1  R    IL  =    ||    [S]"1    [E]    [E]_1  R    \\m 
=    II    [S   -  E]_1    [E]_1   R    || 

—  OO 

■    ■  o  ST«  n  I  jo  'Is  -  «%  tw"1  ^  I 

*   0     "f  ,     jo    I    ([S-Erl»nj    I     II    [«^»IL 

=    (1  +    |    S   -  E    |    +    ...    +    |    (S  -   E)"    |)    ||    [E]"1  R    ||     . 

I  'II  —       I     l  oo 

The  SRC  ensures  that  1  +  ||  S  -  E  |  +  . .  .  +    (S-  E)    is  bounded. 
Q.E.D. 


T 
The  local  discretization  error  I)  =  (d^,  d_  ,  ...,  <ij   is  defined  by 

&o  •  H.(y0;  h)  -  2o    > 

4  -  S  z^  *  fa  lit,,.!,  j^;  h)  -  ^   ,   d  .  1(1)1. 

DEFINITION  3.7.  A  method  satisfying  the  RC  is  said  to  be  quasi- 
consistent   with  the  problem  (2.1)  if 

max    1  ,   1    ~         ,    _ 

..id   M-  0    as   h  -+  0   ■ 
0  $  n  <:  N  '  -n  ' 

and 

max     I   ^    +d+...+d)|+0    as    h  ■*  0. 
0  <  n  $  N  '   —0   —1        -n  .' 

A  method  is   consistent  with  (2.1)    if 

|d^|+|d|    +    ...+    |d_|-»-0         as  h  -►  0 . 

Clearly  consistency  implies  quasi-consistency.   For  a  method  to  be 

quasi-consistent  with  a  problem  it  is  sufficient  that 

(i)   d  =  o(l)    ,   n  =  0(1)N, 
— n 

(ii)  E  d  =  o(h)    ,   n  =  l(l)N 

where  we  employ  the  convention  that  whenever  the  little-oh  or  big-oh  notation 

is  used,  the  implied  limit  or  bound  is  uniform  w.r.t.n.   The  quantity 

Ed  is  the  essential   local  discretization  error. 
— n 

THEOREM  3.8.  A  method  satisfying  the  SRC  is   convergent  for   (2.1)  if 
and  only  if  it  is  quasi- consistent  with   (2.1). 

Proof.      When  R_  =  -  D_,  then  U  =  _Y ;  i.e.,  the  true  solution  is  a  perturbed 
numerical  solution.   (This  unified  treatment  of  discretization  error  and 
roundoff  error  is  due  to  Chartres  and  Stepleman  [l].)   Thus  from  Theorem  3.6 
convergence  is  equivalent  to 

„  ■■   E(cU  +d+...+d  ,  )  +  d    +0    as    h  +  0, 

0  $  n  $  N  '   '—0   —1         -n-1    -n  '  ' 

which  is  equivalent  to  quasi-consistency.   Q.E.D. 
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We  extend  a  theorem  of  Henrici    [3,   p.    71]    stating  necessary  and 

sufficient   conditions   for  convergence. 

THEOREM  3.9.     A  method  satisfying  the  SRC  is  quasi-consistent  with  the 

differential  equation   (2.1a)   if  and  only  if 

(i)  o_(y;    0)    =   Ay, 

(ii)        S   A  =  A, 

(iii)      E  ±{tt    Ay;    0)    =   E   A    f(t,    y) . 

Proof.     Using  the  fact  that  in  the  region  0  $  t  S  1,  0  $  h  <:  h 

continuity  is  equivalent  to  uniform  continuity,  we  have 

(3.10  d^  =  a(70;  0)  -  A  7Q  +  o(l), 

d  =SAy   -Ay  +  0(h) 
-n        n      n 

(3.5)  =  (S  A  -  A)yn  +  0(h),     n  =  l(l)N, 
and 

E  ^  =  E  {v^  +  h  i(tn_l5  v^;  h)  -  yO 

=  E  {^-l  "  ^n  +  h  ^(tn-l'  A  yn-l;  °)}  +  °(h)'    n  =  l(l)N' 
Therefore 

}h   «  4,  =  1  %   -  ^  ♦  jlo  h  1ft,.  A  yj;  0)}  +  0(1) 

t 

=  E  {A(y   -  y  )  +  /n  j,(T,  A  y(  x ) ;  0)  dx}  +  o(l) 
U    n    Z0 
t 

(3.6)  =  E  Q/  n  {-A  f(x,  y(x))  +  i(x,  Ay(x);  0)}dx  +  o(l). 

Assume  (i),  (ii),  and  (iii)  hold.  Then  it  clearly  follows  that 
^0  =  o(l)    , 
d^  =  0(h)    ,   n  =  1(1)N 

n 

,2.  E  d.  =  o(l)    ,   n  =  l(l)N, 
0=1   —J 

which  are  sufficient  for  quasi-consistency.   Now  assume  the  method  is 
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quasi-consistent  with  (2.1a).   Then  d  =  o(l),  n  =  0(l)N,  and  so  from  (3.*0 

and  (3.5)  we  have  a(y •,    0)  -  A  yQ  =  0  and  (S  A  -  A)y  =  0.   Parts  (i)  and  (ii) 

immediately  follow  because  y  is  arbitrary  and  y(t)  is  continuous.   By 

convention  equation  (3-6)  means  that  there  is  some  function  n(h)  tending  to  0 

as  h  ->•  0  such  that 

t 
1  /T*  N  I  E  0X  n  {_A  f(T'  y(T))  +  i(T'  A  y(Th  °)}dT  I  *  n(h)* 

Hence  for  fixed  0  £  t  S  1 

|  E  _/*  {-A  f(x,  y(x))  +  £(t,  A  y(t);  0)}dx  |  <  n( TIT") » 

(J  n  -  q  +  1 

n  =  q,  q  +  1,  •••?  and  so  the  left-hand  side  must  be  zero.   Since  t  is  arbitrary 
and  the  integrand  is  continuous, 

E  {-A  f(t,  y(t))  +  ^(t,  A  y(t);  0)}  =  0 
By  varying  y  ,  y(t)  can  be  made  to  assume  any  value  and  so  (iii)  follows.   Q.E.D. 
For  a  convergent  method  satisfying  the  VSRC  part  (ii)  of  Theorem  3.9 

implies  the  existence  of  a  unique  k  x  s  matrix  M  such  that 

T 
E  =  A  M  . 

The  rows  of  M  are  linearly  independent  left  eigenvalues  of  S  corresponding  to 

the  eigenvalue  1.   Part  (iii)  of  Theorem  3.9  thus  reduces  to 

(3.7)  MT  Mt,  A  y;  0)  =  f(t,  y). 

Remark.      Conditions  (i)-(iii)  of  Theorem  3.9  are  equivalent  to 

convergence  for  methods  which  merely  satisfy  the  RC  if  we  strengthen  the 

smoothness  conditions  on  z_(t;  h)  as  follows: 

(i)   |  _z(t;  h)  -  z_(t;  0)  |  <:  n(h)  where  n(h)  ■>  0  as  h  ■*  0, 

(ii)  z_(t;  0)  is  Riemann  integrable  (cf.  Chartres  and  Stepleman  [l,  p.  U8U]). 

We  can  use  the  fact  that 

I  -10  (S-  E)J| 


IT 

is  bounded  independently  of  n  to  show  that  for  each  problem  (2.1)  there  exist 
positive  constants  c"  and  C"  such  that  for  any  h  S  h  and  R_^  R 

c"  |  |  [E]"1  R  I  1^  <  | |  U  -  U  |  1^  $  C"  | |  R  |  | 
where 

I  I   O   I  I  maX        I    r    T3        I      I        I        v  I  I 

I  I  £  I  I  =  0  S  n  <  N  I  jio  E  ^  I  +  I  ^0  I  +  nil  !  ^n  "  ^n-1  !  ' 
Just  as  before  conditions  (i)  -  (iii)  of  Theorem  3.9  are  necessary  for 
convergence.   On  the  other  hand,  the  additional  restrictions  on  z(t;  h)  make 

it  possible  to  show  that 

N 

I      |  d  -  d    |  =  o(l), 
n=l  '  — n   — n-1  ' 

and  hence  that  conditions  (i)  -  (iii)  are  sufficient  for  convergence. 

Example   I.     Let  us  examine  the  quasi -consistency  conditions  for  a 

straight  multistep  method  satisfying  the  SRC.   We  have 

A  =  e  =  [1,  1,  ...,  1]T 

and  so  part  (ii)  of  Theorem  3.9  requires  that 

(3.8)  -a  -a  -  ...  -a  =1. 

12         q 

Since  the  minimal  polynomial  of  S  is  £  +  ot_  £    +  • • •  +  a  ,  the  SRC  is 

1  q 

T  T 

equivalent  to  the  VSRC,    and  hence  E  =  _e  m     where  m     is  the  left   eigenvalue  of 

T 
S  normalized  so  that  m     e_  =  1.      It  is  easily  verified  that 

m     =p(l)        [l,   1  +  a-.,...,    1  +  a.    +    . .  .   +  a     ,  ] 
—  11  q-1 

where 

p"(0   =  £q  +   (1  +  aj    £q_1  +   ■••   +   (1  +  o.    +   •••  +  a      .). 

1  1  q-1 

Because  i^(t,   u;   h)   =  e_    ^(t,   u;   h),  we  have  from  (3.7)   that 

(3.9)  i|>(t,  A  y;  0)  =  p"(l)  f(t,  y). 

Together  with  part  (i)  of  Theorem  3.9»  (3.8)  and  (3-9)  constitute  the  quasi- 
consistency  conditions  for  straight  multistep  methods.   In  this  case 
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quasi-consistency  is  equivalent  to  consistency  "because  the  local  discretization 

error  e.,  d  and  the  essential  local  discretization  error  p(l)    e  d  are  both 
—In  —  n 

multiples  of  the  same  quantity  d  . 

Example   2.     The  SRC  does  not  exclude  those  weakly  stable  linear 
multistep  methods  based  on  numerical  quadrature.   For  example,  two  steps  of 
Milne's  method, 

u  =u    +  §•  (  \t     +  -t     ,  4f  o)> 

n    n-2   2   3  n   3  n-1   3  n-2 


with  stepsize  h/2  are  equivalent  to  one  step  of  the  method 
1  0 


u 


n-1/  2 


0     1 


n-1 


u 


n-3/2 


+  h 


lf     4f  +if 

6     n        3     n-1/2        6     n-1 


if  4-i^  .    1  - 

6     n-1/2        3     n-1        6     n-3/2 


which  does  satisfy  the  SRC  though  not  the  VSRC.   The  above  pair  of  difference 
equations  is  the  discretization  of  a  system  of  two  identical  differential 
equations,  which  is  not  surprising  because  the  derivation  of  Milne's  method 
involves  two  numerical  integrations  over  each  subinterval. 

Example  S.      The  VSRC  does  not  exclude  linear  multistep  methods  for 
the  special  second  order  differential  equation  y"  =  f(t,  y)  as  long  as  the 
difference  equation  is  written  in  a  form  that  does  not  have  bad  roundoff 
properties.   For  example  the  summed  form  of  Stormer's  method, 

+  h  F 


u  =  u 
n    n-1 


n-1' 


F  =  F  .  +  h  f(t  ,  u  ), 
n    n-1        n   n 


satisfies  the  VSRC. 
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h.      Optimally  stable  methods.   We  wish  to  characterize  those  methods 
that  satisfy  the  most  stringent  definition  of  stability  which  is  compatible 
with  convergence.   Spijker  [6]  has  shown  that  among  convergent  straight  multi- 
step  methods  which  satisfy  the  RC  there  are  some  which  are  more  stable  than 
the  others  and  that  these  optimally  stable  methods  are  those  which  satisfy 
the  strict  root  condition  (which  for  straight  multistep  methods  is  equivalent 
to  the  very  strict  root  condition).   Theorem  U.2  extends  Spijker's  result  by 
considering  a  much  more  general  class  of  methods.   In  the  proof  of  the  theorem 
every  optimally  stable  method  is  shown  to  be  equally  stable  with  a  form  of 
Euler ' s  method . 

LEMMA  k.l.     If  a  method  satisfies  the  VSRC  and  if  a  second  method  is 
at  least  as  stable  as   the  first  method,    then  the  second  method  must  satisfy 
the  VSRC  and  the   E  matrices  of  the  two  methods  must  satisfy   E  E  =  E  . 

Proof.   Call  the  first  method  Method  1  and  the  second  method  Method  2. 
By  Definition  3.2  and  Theorem  3.6  there  is  a  constant  K  independent  of  N  such 
that 

||  U(2)  -  U(2)  ||   .<  K  ||  [E J"1  R  ||  . 

Let  m  be  an  arbitrary  fixed  nonnegative  integer,  and  let  R  =  (iV,  r_  ,  ..., 

T 
r  ,  0 ,  ....  0 )  .   Then 
-m  —       — 

I  "(2)     (2)  I    yr         max    i  „  ,  ,  \ 

u'    -  u      £  K  ..     „     E.  r.+  ...+r  .   +r 
'  -m     -m   '     0  £  n  S  m  '   1  — 0        — n-1    -n 

for  0  <  h  <  min(h  ,  l/(m  +  q  -  l)).   Since  the  left-hand  side  is  continuous 
in  h,  we  can  let  h  ■*  0  to  get 

(k.l)        |  ?nS^"nr  UK*a*    iLlr-.+  M.  +  r  J+  r  |. 
1  n=0  2   -n  '     0  <:  n  <:  m  '   1-0         -n-1   -n  ' 

This  inequality  holds  for  every  nonnegative  integer  m.   Choosing 

r  =  0,  n  =  l(l)m,  shows  that  S  must  be  power-bounded,  which  implies  the  RC 
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is  satisfied.      Suppose,    however,   Method  2  does  not   satisfy  the  SRC.      Then 

there  must  be  an   eigenvalue   £  /  1   of  S     of  modulus   1  and  a   corresponding 

eigenvector  v.      Let 

.n-m 


r     =   £  v        ,        n   =  l(l)m. 

— n  — 


Then  from  (U.l)  we  get 

(m  +  1)  |  v   <:  K 


_n-m    -m 
max    |  £ =_£_  E  v  +  5n-m 

0  £  n  <  m  '     £  -  1    1-   *    - 


2|  E  v  | 

*  k(  i  s  -  i  r  i  2  i}- 

Since  this  inequality  cannot  hold  for  every  m,  we  conclude  that  the  SRC  must 
"be  satisfied.   Therefore  Theorem  3.6  may  he  applied  to  Method  2,  and  so  there 


is  a  constant  c  >  0  such  that 


(k.2)         c  ||  [Ej-1  R  ||   <  ||  UV^  ~  UV 


2)    TT(2 


<  K  |  |  [E  r1   R 


Choose  r  =  (E.  -  I)  v,  n  =  0(l)N,  where  v  is  arbitrary.   Then  (k  .2)   becomes 

— m.  1      —  — 

c  |  N  E  (E  -  I)  v  +  (E  -  I)  v  |  S  K  |  (E   -  I)  v  | . 
Because  both  c  and  K  are  independent  of  N,  E_(E  -I)  v  =  0_,  and  because  v  is 
arbitrary, 

E2  El  =  E2. 
Therefore 

rank  E  =  rank  E  E  £  rank  E   £  s 
showing  that  the  VSRC  is  satisfied.   Q.E.D. 

THEOREM  k.2.     A  convergent  method  is  optimally  stable  if  and  only  if 
it  satisfies   the  VSRC. 

Proof.      Let  a  convergent  method  be  given.  The  global  discretization 
errors  are 

Lq  =  ^y0;  0)  "  A  yo  +  o(1)' 


%  =  s  i(y0;  o)  -  a  yQ  +  o(i). 
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Convergence  implies  that 

£(y;0)  =  A  y, 

S  o(y;  0)  =  A  y 
for  all  y;  hence  S  A  =  A.   This  means  that  the  columns  of  A  are  linearly 
independent  right  eigenvectors  of  S  corresponding  to  the  eigenvalue  1.  There- 
fore there  exists  a  k  x  s  matrix  A  of  rank  s  such  that 

T      T 

Ao s "  V 


Define 


o^y;  h)    =  AQ  y, 


{km3)  so       ■  VAo  Ao)_1  Ao> 

^(t,  u;  h)  =  AQ  f(t,  SQ  u), 

v^(t;  h)    =  AQ  y(t). 
The  method  (4.3)  satisfies  the  VSRC  and  is  convergent  by  Theorem  3.9»   Let  us 
show  that  (4.3)  is  at  least  as  stable  as  the  original  method.   Setting 


6  =  u  -  u  ,  we  have  that 
— n   — n   -n 


^0  =  V 


r=6-S6-h6      ,   n=  l(l)N, 
^i   —n     —n-l     — n-± 


where   6  n    $  L   6  _  I .  Thus 

1  -n-l  '     '  —n-l  ' 

Jo  sno"J  *-rsl±o  +  ih  ST3%  - s  i)-i  - h ij-i>> 

2 
and  since  S     =  S     and  S     S  =  S    , 

I  jo  SV  r-i  I  ■  I  4  +  (so  -  s>  4-i  "  h(  ji  so"J  Ij-i'  I 

«    (1  ♦    I    SQ  -  S    I    +  tn    I    Sq    I    L)    Q     ~,  „    |    ^ 

Therefore  by  Theorem  3.6 
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II  U{0)  -U(0)  ||   (0.  ||  [SJ-1  R  || 

II—         —       I'oo       0  0       —     °° 

.<  CQ  (1  +  I  SQ  -  S  I  +  I  SQ  I  L)  ||  U  -  U  II., 

which  shows  that  (U.3)  is  at  least  as  stable  as  the  original  method.   Having 
established  this  result,  let  us  show  the  necessity  of  the  VSRC  for  optimal 
stability.   Assume  the  original  method  is  optimally  stable.   Then  it  must  be 
at  least  as  stable  as  (U.3),  and  so  by  Lemma  U.l  the  original  method  must 
satisfy  the  VSRC.   Next  we  show  that  the  VSRC  is  sufficient  for  optimal 
stability.   Assume  the  original  method  satisfies  the  VSRC.   Let  any  other 
convergent  method  (Method  2)  be  given  which  is  at  least  as  stable  as  the 
original  method.   We  must  show  that  the  original  method  is  at  least  as  stable 
as  Method  2.   From  Lemma  k.l   it  follows  that  Method  2  satisfies  the  VSRC  and 
that  E  E  =  E  .   Hence  the  null  space  of  E  is  a  subspace  of  the  null  space  of 
E  .   Because  both  methods  are  convergent,  rank  E  =  s  =  rank  E;  so  the  null 
spaces  of  E  and  E  are  identical.   From  E  E  =  E  it  follows  that 

E2(E  E2  -  I)  =  0, 
and  since  E  and  E  have  identical  null  spaces, 


which  simplifies  to 


Therefore 


E(E  E2  -  I)  =  0, 


E  E  =  E. 


n   „    i  ™„v       n 


n-1 


max    |   "   n-j     i      max    ,   "   n-j  ,      )  ji-j 

0  <  n  $  N  '  j  =  0  E    h    '  0  $  n  S  N  '  j=0  E2   -j  U   E2J  j  =  0  E2   ^ 

<  (1    +    I  V  w  M    maX  I   ?   Fn~j  r   I 

$  (1  +  |  E2  -  E  |)0  %   n  g  N  |  .lQ   E2   r_.     [, 

and  so 
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II  H  -  H  I  L  *  C  ||  [E]"1  R  | 1^ 

$  C  (1  +  |  E2  -  E  |)  | |  [E2]_1  R  | | 

^  C  (1  +  |  E0  -  E  | )  c"1  | |  U(2)  -  U(2)  ||  , 

showing  that  the  original  method  is  at  least  as  stable  as  Method  2.   Q.E.D. 


2U 

5.   Order  of  convergence.   In  this  section  it  is  shown  that  the  order 
of  convergence  is  exactly  equal  to  the  order  of  quasi-consistency  if  the  SRC 
is  satisfied.   Then  examples  are  given  of  methods  for  which  the  order  of 
quasi-consistency  is  greater  than  the  order  of  consistency. 

DEFINITION  5.1.   A  method  is  said  to  be  convergent  of  order  p  for  the 
problem  (2.1)  if 

II  U  -  I  I  loo  =  °(hP)  • 
If  a  method  is  convergent  of  order  p  for  all  infinitely  smooth  problems  (2.1), 

it  is  simply  said  to  be  convergent  of  order  p. 

DEFINITION  5-2.   A  method  satisfying  the  RC  is  said  to  be  quasi- 
consistent  of  order  p  with  the  problem  (2.1)  if 

0  .<  n  .<  N  I  ^  I  "  °(h  } 
and 

0  "*i  l«flo*S^"-*V  l  =  °<hP> 

If  this  is  true  for  all  infinitely  smooth  problems  (2.1),  then  the  method  is 
simply  said  to  be  quasi-consistent  of  order  p. 

For  a  method  to  be  quasi-consistent  of  order  p  with  a  problem  it  is 
sufficient  that 

(i)    d^  =  0(hP)    ,   n  =  0(1)N 
(ii)   E  d^  =  0(hP+1)    ,   n  =  1(1)11. 

THEOREM  5.3.     A  method  satisfying  the  SRC  is  convergent  of  order  p  for 
the  problem   (2.1)   if  and  only  if  it  is  quasi- consistent  of  order  p  with   (2.1). 

The  proof  is   omitted  since   it   is    similar  to  that  of  Theorem   3.8. 

Example   I.      If  the  Adams-Bashforth-Moulton  PEC  method  is  written   in 
the   first  way   shown   in  Example  2  of   §2,   we  have 
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E  = 


10      0 
10      0 

10     0 


d     =  h" 
— n 


0 
-5y<3)/12 


+  o(h*), 


and  E  d     =  0(h   ),    shoving  that  the  method  is   quasi-consistent  of  order  three 
~ n 

Written  in  the  second  form,   the  method  is  also   consistent  of  order  three. 
Example  2.     Butcher's  three-stage  procedure  of  order  four  has 

0     0  . 

,        d     =  hk  *°       y     r 

0     1  "* 


E  = 


5f     W   V  ^3> 


+  0(h5) 


and  E  d  =  0(h  ).  There  seems  to  be  no  form  of  this  method  for  which  the 


— n 


order  of  consistency  is  four. 

Example  3.      The  modified  linear  multistep  method  has 


E  = 


1/2 

5/6 

1/2 

1/6" 

"  »ifcW 

0 

0 

0 

0 

k 

0 

1/2 

5/6 

1/2 

1/6 

d  = 
— n 

=  h4 

V>« 

0 

0 

0 

0 

0 

+  0(h5) 


and  E  d     =  0(hP). 
-n 

Example  4.  The  optimally  stable  two-cyclic  method  which  alternates 
Milne's  formula  with  the  third  order  Adams -Moult on  formula  (see  Stetter  [7» 
p.    217])    is   convergent  of  order  four  yet  only  consistent  of  order  three. 
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6.   Asymptotic  behavior  of  the  global  discretization  error.   The 
principal  term  in  the  asymptotic  expansion  of  the  global  error  is  given  in 
Theorem  6.1  for  methods  satisfying  the  SRC.   There  is  some  simplification 
if,  in  fact,  the  VSRC  is  satisfied.   This  is  discussed  after  the  proof  of 
the  theorem. 


THEOREM  6.1.  Assvme   ij>  (t,  u;  h)  is  uniformly  Lipschitz  continuous 

f" 


k 
in  u  and  continuous  in  t  and  h  for  0  s  t  <  1,   u^  R   ,  and  0  £  h  £  h^. 

Assume 


^q  =  hP  y     +  o(hp), 

^  ■  hP  lp(tn)+  hP+1  ^+1    (tn)   +  o(hP+1),  n  =  1(1)N, 

where  <j>   (t)  and  j>     -.(t)  are  continuous  and  E  jL(t)   =   0.      27zen 

u     =  y     +  hP(e(t    )   +   (I   -  S  +  E)"1  6   (t    )   +    (S  -  E)n  c)    +  o(hP) 
— n       jti  —     n  -*p     n  — 

uniformly  for  0  <  n  $  N  where 

c_  =   (I   -  E)   y     -    (I   -  S  +  E)"1  4(0) 
and  the  function  e_  (t)   satisfies 

£(0)    =    E    Y 


e»   =  E  G(t)    (e  +   (I   -  S  +  E)        A    (t))   +  E  A^Ct) 


£■    =  fi  ti^"c  ;    ^_e   +    u    -  b   +  £J~ 

with 

G(t)   =  j^(t,    A  y(t);    0) 
Proof.      Let 


y^  +  h-(e(tn)    +    (I    -  S  +   E)'1  ip(tn; 


u  :      +  hP(e.(t_)    +    (I   -  S  +   E)        4_(0    +   (S  -  E)n  c_). 


Then  we  must  have 


£o  =  Uq  -  i(y0s  h), 


r     =u     -Su     ,    -  h 


,      n  .  i(t      _,  u     _;   h)s  n  =  1(1  )N. 

— n       — n  — n-1  -1-     n-1     — n-1 

It  is  easily  verified  that  _c  has  been  chosen  such  that  r_  =  o(h  ).   By  the 


27 

mean  value  theorem  and  the  Lipschitz  continuity  of  $> 

r     =u    -Su    -   -  h  £ ( t     -, »  y     -, ;  h ) 

-n       — n  -n-1  *■    n-1     lLn-l 

■*  vv  *»-is  h)  (--i  -  ^-i' +  0<h2p+1) 

for  n  ^  1.      It   follows  that 

(6.D       r,  -  -^  +  (^  -  ^)  -  (8  ♦  h  G(Vl))   (Vl  -  ^  ♦  o(hP+1) 

Appropriate  substitutions  yield 

r     =  hP{-4 (t    )    +  e(t    )   +    (I   -  S  +  E)"1   4  (t   )+    (S   -  E)n  c 
-n  -^p     n         —    n  -*p     n  — 


-  S  e(t     .)   -  S(I  -  S  +  E)"1   4   (t     J    -  S(S   -  E)n-1  c}   +  0(hP+1; 
—    n-1  -t?     n-1  — 

Because  E  *   (t)   =  0  and  (I  -  S)    (I  -  S  +  E)        =  I  -  E, 

-4    (t)    +   (I   -  S  +  E)"1   4   (t)    =  S(I    -  S   +  E)"1   4    (t). 

Also  S  e_(t)   =   e_(t)   and  E  c_  =  0_,   and  so   (6.2)  becomes 


r 
— n 


=  hp(e(t    )    -   e(t      .))   +  h*  S(I   -  S  +  E)_x   (4   (t    )    -  4  (t      .))   +  O(h^) 
-    n  n-1  — p     n         -Lp     n-1 

Thus  by  continuity  r     =  o(h    ) .      From  (6.1) 


— n 
— n 


E  r  =  h*"~  E{-4  A,(t  )  +  h   e(t  )  -  h   e(t  ,  ) 
-^p+l  n       —  n       —  n-1 


-  G(t   J  [e(t  -)  +  (I  -  S  +  E)"1  4  (t  .)  +  (S  -  E)n_1  c]}+  o(hP+1)  . 
n-l   —  n-l  ^p  n-l  — 

Making  use  of  the  differential  equation  for  e_(t)  gives 

E  r  =  hP+1  E{h_1  e(t  )  -  h"1  e(t  J  -  e'(t   _)  -  G(t   .)  (S  -  E)n_1  c>  +  o(hP+1 
- n  —  n       —  n-l    —   n-l       n-l  — 


=  -hP+1  E  G(t   .)  (S  -E)11"1  c  +  o(hP+1). 
n-1  — 


Therefore  since 


En         (S   -  E)  <  °° 

n=l    '  ' 


we  must  have 

n-1 


P' 


.Zn   E  r.   +  r  =  o(h^), 

j=0       -j        -n    ' 
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which  by  Theorem  3.6  is  equivalent  to    U  -  U   |   =  o(h  ).  Q.E.D. 

In  the  case  of  a  convergent  method  satisfying  the  VSRC  we  have  that 

T 
E  =  A  M  .   Differentiating  (3-7)  with  respect  to  y  yields 

MT  4  (t,  A  y;  0)  A  =  f  (t,  y), 

and  hence 

E  G(t)  e_(t)  =  E  G(t)  E  e_(t) 

=  A  f  (t,  y(t))  MT  e(t). 

T 
If  we  define  e(t)  =  M  e_(t),  then  the  magnified  error  function  is 

A  e(t)  +  (I  -  S  +  E)"1  6   (t) 


where 


e(0)  =  MT  y 


and 

e'  =  f  (t,  y(t))  e  +  MT[G(t)  (I  -  S  +  E)"1  A  (t)  +  6   +1(t)]. 
Note  that  the  stability  of  the  differential  equation  for  e(t)  depends  only  on 


the  problem  and  not  on  the  method. 

COROLLARY  6.2.  Let  the  hypotheses   of  Theorem   6.1  be  satisfied,   and 

let   0  <  0  £  1  he  given.      Then 

u  =  y  +  hP(e(t  )  +  (I  -  S  +  E)"1  <b  (t  ))  +  o(hP) 
— n   *-n      —  n  -^p  n 

uniformly  for   9  N  $  n  $  N. 

Proof.      For  8  N  $  n  <  N  we  have 

(S  -  E)   c_  !  <:  Bp  for  some  0  <  p  <  1 

_  6N 
$  Bp 

_  -6q   6/h 
<  Bp   H  p    , 

which  goes  to  zero  faster  than  any  power  of  h.   Q.E.D. 

Remark.      Gragg  (see  Stetter  [7,  p.  23^])  has  shown  that  for  p-th  order 

convergent  linear  multistep  methods  satisfying  the  SRC  the  error  posesses  an 
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asymptotic  expansion  of  the  form 


J 


J+l. 


u  =  y(t  )  +  .1     hd  e,(t  )  +  0(h   ) 

n   "   n    j=p     j  n 

uniformly  for  9  H  n  v<  11  as  long  as  f(t,  y)  posesses  J-th  partial  derivatives 
which  are  uniformly  Lipschitz  continuous.   It  is  reasonable  to  expect  this 
result  to  be  valid  for  all  fixed-stepsize  methods  satisfying  the  SRC  as  long 
as  the  increment  function  is  smooth  enough. 

Example   I.     For  the  straight  multistep  method 

^  =  hP+1%Vl(tn)+0(hP+2>- 

By  the  remarks  following  Theorem  6.1  the  magnified  error  function  is 

e_(t)  =  e_  e(t) 
where  e(t)  solves  the  initial  value  problem 

+  , 


e(0)  =  p(1)  1   [en  _  +  (1  +  a.)  e 

q-1  1    q-2 

e'  =  f  (t,  y(t))  e  +  pUT1  ♦  +1(t) 


+  (1+«1+  ...  +a   )  6Q], 


assuming 


u.  =  y(t. )  +  hP  e.  +  0(hP+1),    i  =  0(l)q  -  1. 


Example  2.     For  Milne's  method 


d  =  h' 
n 


r.-i  ..(5) 

1 


2880 


+  0(h°)  , 


j»(t,  u;  0)  = 


^(t,  u;  0)  = 


1/3   2/3 
2/3  1/3 

1/3  2/3 
2/3  1/3 


f(t,  uX) 
f(t,  u2) 


f  (t,  u"1) 

•J 


f  (t,  u  ) 


G(t)  =  f  (t,  y(t)) 


1/3  2/3 
2/3  1/3 
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The  magnified  error  function  is  e(t)  where 


e'  =  f  (t,  y(t)) 

-    y 


This  system  can  be  decoupled: 


e  + 


2880 


1/3  2/3 
2/3  1/3 


1   1 
1  -1 
where 

(.1)--fy(t.x(t))  e1 +  *£$-, 

(e2)'  =  -i  f  (t,  y(t))  e2. 

Example  3.      For  the  modified  linear  multistep  method  equivalent  to 
the  fourth  order  Adams -Moult on  method 


E  = 


1/2   5/6  1/2  1/6 

0    0    0  0 

1/2   5/6  1/2  1/6 

0    0    0  0 
1/U8' 


(I  -  S  +  E) 


-1 


d  =  h 
— n 


k 


0 

-1/U8 

0 


*™  *  *5 


l(t,  u;  0)  = 


-1/80 

0 
17/720 
0 
3/8 
1 
-1/2U 
0 


y(5)  +  h5 
n 


1        1/8      0    -1/2U 
0100 
0    -11/2U    1      1/2U 
0101 
1/128 


1/1+8 

-1/1152 
0 


(k) 

f  (t   ,  y  )  y 
y     n       n       n 

+   0(h6)    , 


W  +      1     1    u     2    ,    1     3,1      Ux 
f(t,   —  u     +u     +  7f  u     +   q-  u    ), 
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jjjt,   u;    0)   = 


3/8 

1 

-1/21+ 

0 


yt.K-»2-K  +  K>r|i  Hi- 


G(t)    = 


3/8 

1 

-1/2U 

0 
The  magnified  error   function  is 


fy(t,y(t))   [I    i    |   I]. 


1 

0 

1 

0 


e(t)    + 


1/U8 

0 

-1/1+8 

0 


y(u)(t) 


where 


and 


e(0)   =  I     1  +  1     2  +  I  y3  +  1     J» 
e^o;   "  2  Yl|  +   6  Yl+        2  Yl+        6  YU 


e«   =  f  (t,  y(t))    e  +-^-y(^(t)   +  r4  f  (t,  y(t))  yU)(t) 


(M, 


180 


1+8  V 


This  can  be  rewritten 


(e  +-^y(1°(t))'   =  f   (t,   y(t))    (e  +J^y{k)(t))   +  ^y(5)(t 


1+8 


720 


hence  e(t)  +  rs-  y .   (t)  satisfies  the  differential  equation  for  the  magnified 


error  function  of  the  fourth  order  Adams -Moulton  method. 
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